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A FINITE -DIFFERENCE ANALYSIS OF THE 
NOZZLE STARTING PROCESS IN 
AN EXPANSION TUNNEL 

K. James Weilmuenster 
Langley Research Center 

SUMMARY 

A finite-difference solution of the time-dependent, quasi-one-dimensional equations 
for a bigas mixture, including real gas effects, for both viscous and inviscid flow has been 
applied to the starting process in an expansion tunnel nozzle. 

An analysis of the numerical technique has shown that using a Richtmyer differencing 
solution of the inviscid quasi-one-dimensional equations along with a smoothing function 
avoids numerical instabilities unrelated to classical stability criteria and that the stability 
restriction on step size for the quasi-one-dimensional equations is identical to that for 
the one-dimensional equations. 

The solutions obtained for the starting process in a conical nozzle 1.732 m in length 
indicate that expansion tunnel operation without a ternery diaphragm is restricted to those 
conditions for which the test flow is not disrupted by an upstream -facing shock and that the 
use of a ternery diaphragm in the system , combined with a sufficiently low initial nozzle 
pressure, will allow the nozzle to be started in an acceptable manner. 

INTRODUCTION 

The development and analysis of the expansion tube, a device used to produce high- 
stagnation enthalpy flows, have been thoroughly discussed in several papers (refs. 1 to 4). 
The expansion tunnel is a modification of the expansion tube which utilizes the basic expan- 
sion tube configuration with a ternery diaphragm and a conical nozzle added at the down- 
stream end of the acceleration chamber as illustrated in figure 1. 

Trimpi and Callis (ref. 5) have made a perfect gas analysis of the expansion tunnel. 
Their calculations determined, among other things, the test time lost as a result of the 
starting and stopping processes in the nozzle, as well as those conditions which must be 
satisfied if there is to be a perfect nozzle start (e.g., that no upstream compression waves 
be produced by the passage of the primary shock through the nozzle (see ref. 6)). 


Friesen and Moore (ref. 7) measured velocity profiles and pitot pressure in the exit 
plane of an expansion tunnel nozzle. They investigated the two possible modes of expan- 
sion tunnel operation. First, the facility can be operated without a ternery diaphragm, 
which results in the gas and initial pressure in the nozzle and dump tank being the same 
as that found in the acceleration chamber. Second, the ternery diaphragm can be used, 
which allows the gas and initial pressure in the nozzle and dump tank to be varied. Using 
the shock in the acceleration chamber to burst the ternery diaphragm is a good approxima- 
tion of the assumption of an instantaneous diaphragm opening used in the theoretical anal- 
ysis of the expansion tunnel. However, experiments have shown that allowing the incident 
shock in the acceleration chamber to burst the ternery diaphragm may degrade the quality 
of the flow entering the nozzle. Therefore, Friesen and Moore employed an electromag- 
netically opened diaphragm whose construction and effect on the flow entering the nozzle 
are detailed in reference 8. Essentially, they found that it was necessary to use a ternery 
diaphragm and to evacuate the nozzle before reasonably uniform nozzle exit conditions 
could be found. 

The results of reference 5 are limited to those initial conditions for which no 
upstream-facing shock is generated in the nozzle. In order to extend the analysis of 
expansion tunnel operation to conditions which generate shocks in the nozzle and to make 
parametric studies of the effect of these shock waves on the test flow quality, it has been 
necessary to develop a shock- capturing, finite -difference method for solving the funda- 
mental equations of quasi-one-dimensional flow. 

A finite-difference method for solving the time-dependent, quasi-one-dimensional 
equations for a real, bigas flow is investigated for both viscous and inviscid flow. The 
method deemed to be most appropriate has been used to analyze the starting process in an 
expansion tunnel. Specifically, the effect of initial nozzle pressure and the formation of 
internal shocks on the starting process and the available test time are investigated. 


SYMBOLS 

A nozzle cross-sectional area, m2 

A area ratio, A(x)/A(x = 0) 

a sound speed, m/s 


B,C,F,J,'^ 
K,L,Q J 


coefficient matrices 


specific heat at constant pressure, m2/s2-K 
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D 

E 

e 

f,g,s 

i 

k 

Ms 

Npr 

NRe 

P 

Rg 

T 

t 

At 

U 

u 

W 


specific heat at constant volume, m^/s^-K 
nozzle diameter, m 

.cal energy, p(e . kg/™-s2 

internal energy, m^/s^ 
vector-valued functions 

thermal conductivity, N/s-K 

incident shock Mach number in region (I) 

molecular weight of species s, g mole 

Prandtl number 

Reynolds number 

pressure, N/m^ 

gas constant of species s, m2/s2-K 
temperature, K 
time, s 

time increment, s 
free- stream velocity, m/s 
velocity, m/s 

vector of dependent variables 


w 


mass fraction 


X distance, m 

Ax spatial increment, M 

y ratio of specific heats 

6 nozzle half -angle, deg 

€ error in solution as function of time and space 

X mesh ratio. At/ Ax 

/i viscosity, N-s/m^ 

u eigenvalue of amplification matrix 

p density, kg/m^ 

Subscripts: 

(T) region (see fig. 1) 

j part of a sequence 

mix property of gas mixture 

0 initial or reference value 

s particular species 

Superscript: 

^ dimensional quantity 

Abbreviation: 

SSS secondary starting shock 
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STATEMENT OF PROBLEM 


The manner in which the starting flow in an expansion tunnel nozzle develops, for a 
given set of test conditions, has a direct effect upon the operation of the facility. Specifi- 
cally, the experimental results of Friesen and Moore (ref. 7) indicated that an electromag- 
netically opened ternery diaphragm had to be included in the system before a uniform test 
flow could be established at the nozzle exit. The addition of this hardware increases the 
system maintenance time between timnel runs as well as adding further complexity to the 
facility operation. Also, much of the facility turnaround time is spent evacuating the var- 
ious sections of the system of which the nozzle dump tank is by far the largest. Obviously, 
higher allowable initial nozzle pressures will require shorter pumping times. Since 
Friesen and Moore maintained the initial nozzle pressure at approximately 10"^ torr 
(1 torr = 133.3 N/m^), the effect of initial nozzle pressure on the nozzle starting process 
is unknown. 

In this report, a numerical method is developed to determine what usable test flows 
can be generated, for the given nozzle geometry and initial entrance conditions, without 
the use of a ternery diaphragm, and what is the maximum allowable initial nozzle pressure 
when operating with a ternery diaphragm in the system and the penalty for exceeding this 
value. 

The paper deals with the starting process in a conical nozzle having an entrance 
diameter of 0.076 m and a length of 1.609 m, whose entrance conditions correspond to 
those found in the acceleration chamber of an expansion tube (i.e., a low-Mach-number 
flow followed by a high-Mach-number flow). These dimensions conform to those of a 
nozzle presently installed in the Langley 6 -inch expansion tube. In the operation of an 
expansion tunnel, three different gases come into play: the test gas, initially in the inter- 
mediate chamber, the acceleration gas, initially in the acceleration chamber, and the gas 
initially in the nozzle, all of which may be the same or different. In this study, only two 
different gases will be used at any one time. Thus, it will be necessary to include the 
dynamics and thermodynamics of a bigas flow in the analysis. 

A wave diagram for the flow that might be expected for an imperfectly started noz- 
zle is indicated in figure 2. The gas initially in the acceleration chamber, which is pro- 
cessed by the primary shock, and the gas initially in the nozzle may or may not be identical 
or have identical pressures, depending on whether or not there is a diaphragm separating 
the nozzle from the rest of the system. For the computations in this report, it is assumed 
that if a ternery diaphragm does exist in the system, the diaphragm opens instantaneously 
upon contact with the incident shock. In the nozzle, the strength of the primary shock is 
decreased by the increasing area ratio while an interface separates the shocked 
acceleration- chamber gas from the shocked nozzle gas and an unsteady expansion fan. 
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This region of unsteady flow is terminated by an upstream-facii^ secondary starting 
shock. The secondary starting shock is followed by the steady expansion of the shocked 
acceleration gas. As the interface separating the shocked acceleration gas and the test 
gas traverses the nozzle, it is accelerated and is preceded by a family of downstream- 
facing characteristics and followed by a family of upstream-facing characteristics. When 
the last of these has passed the nozzle exit, the nozzle is fully started. Test time in the 
nozzle is terminated by the arrival of the trailii^ edge of the expansion tube expansion fan 
at the nozzle exit. 


METHOD OF SOLUTION 

Several investigations have been made of the starting process in a hypersonic noz- 
zle, most of which were performed in reflected shock tunnels. Perhaps the best known of 
these studies was published by Smith (ref. 9), who used the method of characteristics and 
a perfect gas to obtain solutions. More recently. Stalker and Mudford (ref. 10) made an 
experimental and simplified analytical analysis of the starting process in a nonreflected 
shock tunnel. 

The analysis of hypersonic flows in nozzles by using finite -difference solutions was 
first investigated by Crocco (ref. 11). Crocco utilized the quasi-one-dimensional equa- 
tions in his work which concentrated on the stability of such a system rather than on the 
solution of any particular problem. Laval (ref. 12) used a two-step, second-order- 
accurate differencing scheme that employed an explicit artificial viscosity along with the 
quasi-one-dimensional equations to solve the reflected shock tunnel starting problem and 
obtained results which were in good agreement with those of Smith (ref. 9). Laval also 
used a perfect gas. Benison and Rubin (ref. 13) used a variation of the Richtmyer method 
to solve the quasi-one-dimensional equations for a perfect gas including the effects of 
viscosity and heat conduction. They, however, limited their calculations to flows contain- 
ing only weak shocks. 

Since the end result of this study is to provide time histories of nonlinear distur- 
bances as they traverse a nozzle, the finite-difference scheme used should be second- 
order accurate in both time and space to insure good space -time resolution. A survey of 
available methods led to the choice of the Richtmyer differencing scheme as developed by 
Thommen (ref. 14) for use in this analysis. This method closely resembles the Benison- 
Rubin method yet requires fewer numerical computations. Also, there is some ambiguity 
associated with the specification of the upstream boundary conditions in the Laval method. 
In addition, the effects of imperfect gases and bicomponent flow, which would unduly com- 
plicate a characteristics solution, are included in the present analysis. 
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Briefly, the Richtmyer method is as follows. Given the partial differential equation 
of the form 


Wt = -f{W)x + s(Wx), 


( 1 ) 


the Richtmyer method takes on the form 


W(t + At,x) = W(t,x) - —I 

Ax 




+ At 


S(W), 


( 2 ) 


where 


W 


t + y,x ± = |[w(t,x) + W(t,x ± Ax) T f(t,x) - f(t,x ± Ax) 


At 


N 


x±Ax 


XX 


xy 


(3) 


This method is the result of a Taylor series expansion of W about the point x 
and gains its second-order accuracy through the use of an intermediate step. 

Details of the Richtmyer method are developed as follows. Expand W(t,x) about 
the points wft + ^,x ± ^ such that 


^ " f’" ^ " f r f " f " 


(4) 




= ^ - f j " f - f - 0(At)2 


(5) 


and since W|. = -f^ + Sx 


^ " f " f = " f - f " f r f " fi " 
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(7) 


W(t + f,x - = W[t,x -fj-f fxp - #) + Sx(tpc - + 0(At)2 


Ax I ^ At 

2 / 2 


Then, if W is expanded about the point W(t+^yXj, 

a 


W(t + At,x) = wft + |^,xj + I* Wtft + |^,xj + Wttft + |i,xj + 0(At)3 


( 8 ) 


W(t,x) = W t + |^,x - Wtft + |^,xj + ^ Wttft + |^,xj + 0(At)3 


(9) 


Subtracting equation (9) from equation (8) and rearranging gives 


W(t + Atpc) = W(t,x) + AtWt^t + y,x| + 0(At)^ 


( 10 ) 


or with the substitution Wt = -fx(W) + S(W)jqj 


W(t + At,x) = W(t,x) - Atfxft + ^,xj + AtSxft + |t,xj + 0(At)2 (11) 


The spatial derivative f^^t + is evaluated by using the values of W at 

|t + At,x ± To maintain the formal second-order accuracy of equation (11), second- 

order-accurate central differences are used for the spatial derivatives in equations (4), 
(5), and (11); that is. 


and 


fx t,x ± ^ = 


Ax \ _ Tf(t,x) ± f(t,x ± Ax) 


2 (Ax/2) 


Ixfl + 


At Ax\ ,/. At Ax 
2(Ax/2) 
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The terms 



in equations (6) and (7) are the averages 


|[w(t,x) + W(t,x ± Ax) 


Sx in equation (11) is evaluated at (t^) rather than at ^t + improve the computa- 

tional efficiency of the scheme. This change has been used by Thommen and several other 
authors without any apparent effect on computed results. 

The quasi-one-dimenslonal equations, which describe the flow of a multicomponent 
gas including the effects of viscosity and heat conduction, can be written in the followir^ 
manner in the vector form 

Wt = -fx - Agx + s(Wx)^ (12) 

where 



pA 


puA 


0 


0 


pWgA 


^ ^ ^ 

puAwg 

II 

0 


0 

w = 

pGa 

f = 

— 2- 
pu A 

1 AuxA 

g = 

p 


pEA 

i 


uPa(e + 


pu^juA + AkT- 


0 


and E = e + ^ u2, w„ = A = A(x), ^ = m(T), k = k(T), and p = p(e,p). The inviscid 

2 p 

forms of the equations are obtained by eliminating the S matrix from equation (12). 
These equations are nondimensionalized by using the following definitions: 


u 

D-1- 

T 

T - 


X 

D - P 

U — — 
Uq 

^■Po 

T 

' 'lo 

A. “ ■ 

^o 

^ Po^oTo 

h 

Ws = — 
p 


IX - — 

Mo 

k=l 

ko 

A = ^ 

X 2 

*o 


where Xq 

is taken to be the square root of the minimum 

value of 

A. The resulting 


equations are 


9 



4 


(14) 


(pA\ = -(puA)jj 

(pAWs)^ = -(puAws)j 


(puA)t = -(pu2Al-^^P: 


TT 2 X 3N-« 
Uq^ '^^Re,o 


(mAUx) 


(PEA)^ = 


pAuf^^P 

\P py 


3N, 


Re,o 




Since equations (14) are set up to calculate bigas flows, the transport properties for a 
bigas mixture must be determined. Equations and methods for this procedure are pre- 
sented in appendix A. The temperature and pressure of the bigas mixture are determined 
by an iterative procedure using the equation for the internal energy of a gas mixture. 
Helium is taken to be a calorically perfect gas and the real air properties are obtained 
from the subroutine of Tannehill and Mohling (ref. 15). 


Stability Analysis 

A classic Von Neumann stability analysis (ref. 16) of equations (14) in a linearized 
form will be made. 

The linearized equations (14) can be written as 


Wt = -(B + AF + Q)Wx + CWxx 


(15) 


where B = C = , Q = and F = Benison and Rubin worked out these 

8W aWx 8W 8W 

Jacobians allowing for the temperature dependence of p and k while assuming the 
working fluid was a perfect gas. The final form of their amplification matrix was of such 
complexity that closed-form solutions for the eigenvalues could not be found. 

In an effort to find closed-form solutions for the stability limits, the followii^ anal- 
ysis was made. Except for the form of the equations used, this analysis parallels that of 
Thommen (ref. 14), who found the stability criteria for the one -dimensional viscous equa- 
tions by usii^ the Richtmyer method. Assume that in the vicinity about a point in the flow, 
the viscosity, thermal conductivity, and ratio of specific heats, although functions of tem- 
perature, are constant. As discussed by Richtmyer and Morton (ref. 17), the Von Neumann 
stability condition for the nonconservative form of the governing equations is the same as 
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that for the conservative form. If the linearized form of equations (14) is written in the 
nonconservative form, the Jacobians will be greatly simplified when compared with those 
used by Thommen, who considered only the conservative form of the equations. This sim- 
plification makes the determination of the eigenvalues of the amplification matrix much 
easier. 

If the definitions p = (y - l)pe and 2 ? - y(y - l)e are used, and the momentum 
equation is redefined as (puA)t = -(pu^A + pA) + pAx + 


3^Re,o 


(pAuxL> equations (14) 


can be written in their nonconservative form as 


- 





“ 

- 

PA 


-u 

0 

-P 

0 

pA 

wA 


0 

-u 

0 

0 

wA 

uA 


®v,o'^o^^ 

0 

-u 

(7o “ ^)^V,o'^ 

uA 



eA 

t 

0 

0 

y 

-u 

eA 


PU 


0 


^v,o'^o^^ 

Ajj + 

r 

_ ua^ 


y 



0 

0 


0 ^ 


0 
0 

4 p/p 


3 N- 


Re,o 

0 


0 

0 

0 

y(p/p) 


NRe,oNpr,o 


pA 

wA 

uA 

eA 


(16) 


-ixx 


or 


Wt = JWx + KAx + LWxx 


(17) 


where the coefficient matrices J, K, and L are taken to be constant for the stability 
analysis. Now, by using forward time and central space differencing, equation (17) can 
be written as 

W(t + Atpc) = W(t,x) + j[w(t,x + Ax) - W(tpc - Ax)] + k[a(x + Ax) - A(x - Ax)] 


At 


2 Ax-^ 


L[w(t,x + Ax) - 2W(t,x) - W(tp£ - Ax)] 


(18) 
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Following reference 18, assume W(t,x) is an exact solution of equation (18) and e(tpc) 
is the error in the solution at x and t. Then, substituting W(t,x) = W(tpc) + e(t,x) into 
equation (18) gives 

W(t + At,x) + e(t + At,x) = W(t,x) + e(t,x) + j[w(t,x + Ax) + e(t,x + Ax) - W(t,x - Ax) 


- e(t,x - Ax)J + K^A(x + Ax) - A(x - Ax)j + — ^^[w(t,x + Ax) 
+ e(t,x + Ax) - 2W(t,x) - 2e(t,x) + W(t,x - Ax) - e(t,x - Ax)j (19) 


When equation (18) is subtracted from equation (19), an expression for the error at 
(t + At,x) is found, and the expression 


e(t + At,x) = e(t,x) + 


At 

2 Ax 


e(t,x + Ax) - e(t,x - Ax) + L e(t,x + Ax) 

2 Ax^ •- 


- 2e(t,x) + e(t,x - Ax) 


is independent of the expression KA^. Therefore, only the stability of the expression 


Wt = JWx + LWxx 


( 20 ) 


will be investigated. 

For the inviscid case, equation (20) reduces to = JWx, the stability criterion of 
which is established by the well-known CFL (Courant-Friedrich-Lewy) condition 
At/ Ax = l/(u + a). 

Following Von Neumann, equation (20) is differenced according to the Richtmyer 
method. However, for the case of constant coefficients, the Richtmyer method reduces 
to the original Lax-Wendroff method; that is, 

W(t + At,x) = W(t,x) + jrw(tpc + Ax) - W(tpc - Ax)l + j2[w(t,x + Ax) 

2 Ax L J 2 Ax^ 


- 2W(t,x) + W(t,x - Ax) 


JL At^ 
2Ax3 



+ Ax) - Wxx(t,x 



+ L At Wxx(l5x) 


where higher order terms in L have been neglected. 


( 21 ) 
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If a Fourier component of the solution e®*e^^ is introduced in equation (21), and 

0 !(t+At) i)3x 

the amplification matrix is defined as , the following expression is obtained: 

„Q!t i/3x 
e e ^ 


gOfAt _ j ^ cos(^ Ax) - ij + Ji sin /3 Ax 

where X = At/Ax. A necessary condition for the stability of the error e®^^ is 

jeOtAtj g Following Thommen, the largest values of the amplification matrix occur 
when ^ Ax = IT which reduces equation (22) to 


X + 


2LX^ 


Ax 


cos(/3 Ax) - 1 


( 22 ) 


ofAt 

e 


x2j2 



(23) 


The determination of the eigenvalues of equation (23) is a long, tedious exercise in 
algebratic manipulation and will not be presented here. The results are identical to those 
found by Thommen who used the conservative form of the one-dimensional equations. 
Thommen^s analysis did not include the species concentration equation. This, however, 
only resulted in a repeat of one of the eigenvalues found when the species -concentration 
equation was excluded from the analysis. The eigenvalues for the case where heat con- 
duction is neglected are 


1 - 2X^u2 
1 - 2X^u2 


1-2 




x2(u2 + a2 


4Xp. 


PNr 


e,o 


Ax 




16x2/i2 


1/2 


^2 t , t 2 * 5 

P %e,o 


1 - 



4Xp 

^ P%e,o ^ 


(4x4u2a2 + ^ ^ 


1/2 


(24) 


Attempts were made to determine the eigenvalues for the case where heat conduc- 
tion was included. The species-concentration equation was eliminated so that the charac- 
teristic equation could be simplified. The third-order equation which must be solved for 
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the eigenvalues is quite involved in this case. Therefore, the eigenvalues of the amplifi- 
cation matrix were found by using a symbolic manipulator computer code. The resultir^ 
answers are presented in appendix B. One is real, the other two are complex. The max- 
imum eigenvalue will most likely come from one of the complex forms. Since the condi- 
tion for stability requires the magnitude of the eigenvalue to be less than one, it would be 
necessary to evaluate the roots of ^ to find the appropriate value of X. It would 

be preferable if such a calculation could be avoided by considering (1) the complexity of 
Vj , (2) the number of times ^ would have to be evaluated, and (3) that an iterative 

procedure would have to be used to determine X. Thommen (ref. 14) looked at the sta- 
bility of the equations in a quiescent region when heat conduction was included and found 
the stability criterion to be too restrictive. He obtained good results in subsequent cal- 
culations by basing the stability on the eigenvalues (eq. 24). Thus, the stability of the 
viscous equations will be based on the eigenvalues (eq. (24)), 


Solution Stability 

To test the computational scheme, a simplified test case was run for a nozzle having 
p^^^ = p^g^ = 16.67 N/m2 in helium, a primary shock Mach number of 5, and no ternery 

diaphragm in the system. 


Equations (14) were used in their inviscid form which resulted in a numerical insta- 
bility just upstream of the secondary starting shock. Figure 3 shows the evolvement, in 
time, of the axial distributions of velocity and log^Q of the density for this test case. 

The continuous upper and lower portions of the u/uq and logjQ curves, respectively, 

^o 

represent the steady expansion solution for the nozzle. The position of the primary shock 
is indicated on the velocity curve for the first time interval shown. Inspection of the first 
curve shows no evidence of a secondary starting shock (SSS). The second curve shows a 


steepening of the u/uq and 



gradients near the steady expansion curve, while 


the third curve indicates a definite shock at the same point accompanied by a slight over- 
shoot in the velocity and density. As the SSS traverses the nozzle, oscillations associated 
with the numerical solution appear at the upstream side of the shock. These ultimately 
cause a negative temperature to appear, leading to a rapid degeneration of the solution. 


The instability in the solution did not seem to be related to the step size used, since 
the oscillations at the shock had appeared when it was first formed near the nozzle 
entrance and did not grow as the shock moved downstream. Only when the temperature at 
the shock, because of the highly expanded nature of the flow, became small when compared 
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with the oscillations in the temperature did the negative temperature and subsequent insta- 
bility occur. Since all second-order -accurate and higher order accurate difference 
schemes have oscillations at a shock unless sufficient damping is provided (e.g., artifi- 
cial viscous terms used by Laval and others), the full viscous form of equations (14) was 
used to see if the real viscous effects would be sufficient to damp the oscillations at the 
shock. This seemed to be a reasonable course to follow since the free-stream unit 
Reynolds numbers were small enough to allow cell Reynolds numbers in the range of 2 
to 10. 


Sufficient smoothing of the solution was accomplished by using the viscous form of 
equations (14) and a spatial increment small enough to permit resolution of the secondary 
starting shock. In this manner, satisfactory results were obtained, although some accu- 
racy had to be sacrificed as a result of the smearing of shocks and interfaces by the vis- 
cous effects. An example of these results is shown in figure 4. As the flow progresses 
through the nozzle, the shock is smeared out and the overshoot disappears. Note, how- 
ever, that the oscillations at the primary shock never disappear. The cause of this effect 
is readily apparent when the nonconservative equations (eqs. (17)) are examined. The 
coefficients of the Wxx terms are made up in part by a term M/pNRe,o> which is, in 

effect, an inverse unit Reynolds number based on the local free-stream density and vis- 
cosity and on the reference velocity. The value of this Reynolds number varies widely 
along the nozzle as the solution evolves in time. Small values of the inverse Reynolds 
number force the equations to the hyperbolic limit, for which oscillations at discontinu- 
ities are assured when second-order-accurate difference schemes are used. Large val- 
ues, on the other hand, maintain the parabolic character of the equations and smear out 
discontinuities if the spatial increment is sufficiently small to provide good resolution. 
Thus, for the case shown in figure 4, equations (14) are parabolic-dominant in the vicinity 
of the SSS and hyperbolic-dominant near the primary shock. 


Although the use of the viscous form of the equations to stabilize the SSS worked 
well in this case, the very nature of expansion tube flow suggested that such a procedure 
would be incapable of handling the range of conditions required for this study. In refer- 
ence to the following sketch, the local Reynolds number in region (2^ is much less (an 
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order of magnitude or more) than the local Reynolds number in region Thus there 
would be a large jump in the inverse Reynolds number across the interface. Initial cal- 
culations showed that if the SSS followed path A, it would remain stable, whereas, if the 
shock followed path B and moved across the interface, the computation would destabilize 
at the SSS. This behavior is consistent with the previous observation of the effect of the 
local inverse Reynolds number (i.e., cell Reynolds number) and its relation to the stability 
of the SSS. 

Thus, to generate a general algorithm for this problem, it was necessary to intro- 
duce a stabilizer for the SSS. One option would have been to implement the stabilizers 
used by Richtmyer (ref. 16) or Laval (ref. 12), or to decrease the cell Reynolds number 
in the expanded gas of region by a corresponding reduction in grid size. Use of the 
inviscid form of the equations along with the simple smoothing method proposed by 
Lapidus (ref. 19) appeared to be a more attractive option. 

The Lapidus method involves the solution of a diffusion- like equation involving the 
conserved flow quantities in the updated time plane. Given the equation 


Wt = -f(W) 

with a solution W(t + At,x) in the updated time plane, the smoothing is accomplished by 
solving the equation 


W(t + At,x)j^g.^ = W(t + At,x) + ri( |u(t + At,x + Ax) - u(t + At,x) 


W(t + At,x + Ax) 


- W(t + At,x) 


- W(t + At,x - Ax) j 


u(t + At,x) - u(t + At,x - Ax) 




W(t + At,x) 


where W is the vector defined in equation (13). Simple, one -dimensional shock compu- 
tations using the smoothing function indicated that rj = 0.5 provided sufficient smoothing 
of preceding postshock oscillations without causing excessive diffusion of the shock front. 
Results from inviscid test case computations which included the Lapidus smoothing are 
shown in figure 5. When compared with the viscous solution (fig. 4), there is a persistent 
small overshoot at the SSS, but the shock resolution is much better than that found in the 
viscous case. Subsequent calculations showed this method would maintain the stability of 
the SSS over a wide range of conditions and it has been employed in the following compu- 
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tations. Use of the Lapidus smoothii^ does not affect the stability criterion for the 
numerical method which was taken to be the standard CFL condition for the inviscid form 
of equations (14). 


RESULTS 

Nonevacuated Nozzle 

Inclusion of a ternery diaphragm in the expansion tunnel system is based on the 
assumption that to have acceptable nozzle starts, the nozzle would have to be at a much 
lower initial pressure than the acceleration chamber. Since operation of the facility 
would obviously be easier if this diaphragm could be eliminated, the first cases were run 
to determine if acceptable nozzle starts could be made when the nozzle pressure was the 
same as the acceleration chamber pressure. 

Figure 5 shows data for a case where an initial shock Mach number was 5 at the 
nozzle entrance and where the arrival of the interface separating the acceleration gas and 
test gas was not taken into account in the computation. In determining the conditions to be 
used for region (^, the pressure and velocity in regions and must be matched. 

The matching is done by finding simultaneous solutions for the Rankine-Hugoniot relations 
across the shock in the acceleration chamber and by finding the relationships governing 
the unsteady expansion of the flow from region (^. Typical expansion tube conditions of 
50 torr in air in region (T) and a shock Mach number of 7.93 in the intermediate chamber 
were used to obtain the conditions for region (^. Matching the conditions required a 
shock of = 7.95 moving into region which had an initial pressure of 0.886 

torr in helium or a shock of Mg = 16.59 moving into air at 0.360 torr. The initial 
conditions and reference values used in these computations are in table 1. The time 
interval between shock and interface arrival at the nozzle entrance was taken from exper- 
imental results to be 50 /ns and 20 /ns for the helium and air acceleration gases, 
respectively. 

Figure 6 shows some partial results for a 10° nozzle and helium acceleration gas. 
Pertinent facets of the flow can be determined from such plots and transfered to x-t wave 
diagrams for more convenient comparison and analysis. Cases were also run for an 8° 
nozzle with a helium acceleration gas and a 10° nozzle with an air acceleration gas. 
Results for these three cases are shown in figures 7, 8, and 9. In all three cases, the SSS 
gains strei^th quickly due to the large back pressure in the nozzle. After the SSS has 
passed through the interface, it is swept downstream at a more rapid rate because the 
free-stream Mach number of the test gas is greater than that of the acceleration gas. At 
approximately 0.6 m downstream of the diaphragm, the SSS and the trailing edge of the u-a 
family of starting waves intersect. After this intersection, the SSS moves downstream at 
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a decreasing rate, thus affecting a portion of the test gas. In each case, substantial losses 
in available test time are indicated at the larger values of A. 

This conclusion is in agreement with the experimental observations of Freisen and 
Moore (ref. 7), who found the test gas flow quality for nonternery diaphragm conditions to 
be very poor when compared with the flow observed when a ternery diaphragm, which 
results in a lower nozzle pressure, was used. 


Evacuated Nozzle 


Inclusion of a ternery diaphragm in the system allows the initial nozzle pressure to 
be varied. A reduction in the value of p^ should lead to a reduced influence of the SSS 

on the test flow time. The three previous cases were recomputed by using pressure 
ratios across the ternery diaphragm of 0.01. The resulting x-t 

diagrams (figs. 10 to 15) show that in each case, subsequent reductions in P^^^P^^ 


led to reductions in test-time losses. 


The effect of P^^y^P^^^ on the available test time is summarized in figures 16 

to 18 for each of the three cases. For both the 8° and 10° cases with helium as an accel- 
eration gas, approximately 0.01 is needed to eliminate the effects of the 

SSS on the test gas flow, whereas a of only 0.1 is needed when air is used as 

the acceleration gas. 


The time of 200 jis between the arrival of the interface and the expansion fan at the 
nozzle entrance is based on experimental observations. Theoretically, this time could be 
increased by lengthening the acceleration chamber, thus decreasing the effect of the SSS 
on the available test time. However, viscous effects, as detailed in reference 20, degrade 
the quality of the flow in region after it has traveled a sufficient distance downstream 
of the secondary diaphragm. The 200 jis of test time used in this case is near the maxi- 
mum amount of good quality test flow that might be expected for these flow conditions. 


Practical test setups require a balance between test time, test Mach number, and 
model size. Also, in this case, one of the factors in testing frequency is the ratio 

For some running conditions, P^^ very small; the facility configuration 

is such that the volume of the nozzle and dump tank is very large, which for small values 
of p^ leads to long nozzle evacuation times. To evaluate these effects for the data 

which have been presented, a plot of test time against exit Mach number is shown in fig- 
ure 19 and a plot of test time against model diameter, based on the nozzle entrance diam- 
eter, is shown in figure 20. The test flow Mach number can be doubled and a model size 
four times the nozzle entrance diameter is available in a system which excludes a ternery 
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diaphragm without suffering SSS-related test-time losses. Further increases in exit Mach 
number and/or model size will result in drastic reductions in available test time. This is 
also the case when a ternery diaphragm is used with helium acceleration gas and when 
p^yp^^^ = 0.1 is used for either nozzle configuration. Test flows are not free of inter- 
ference from the SSS until 

p /p/-^ = 0.01 is used with a helium acceleration gas. 

®/ @ 

CONCLUDING REMARKS 

A finite-difference solution of the time- dependent, quasi -one -dimensional equations 
for a bigas mixture, including real gas effects, for both viscous and inviscid flow has been 
applied to the starting process in an expansion tunnel nozzle. 

An analysis of the numerical technique has shown that using a Richtmyer differencing 
solution of the inviscid quasi-one-dimensional equations along with a smoothing function 
avoids numerical instabilities unrelated to classical stability criteria and that the stability 
restriction on step size for the quasi-one-dimensional equations is identical to that for 
the one-dimensional equations. 

The solutions obtained for the starting process in 8° and 10° conical nozzles having 
an entrance diameter of 0.076 m and a length of 1.609 m show that expansion tunnel oper- 
ation without a ternery diaphragm is limited to moderate Mach number test flows (Mach 
number of 12 to 13) and small area ratios (for example, small model diameters). The 
solutions also indicate that the use of a ternery diaphragm and ternery diaphragm pres- 
sure ratios of 0.1 and 0.01 for air and helium acceleration gases, respectively, allows the 
full potential of the nozzle to be realized (that is, expansion to a Mach number of about 18 
in the nozzles with an area ratio of 70). 

Langley Research Center 

National Aeronautics and Space Administration 
Hampton, Va. 23665 
November 14, 1975 


(N)/® 


= 0.1 is used with air acceleration gas or until 
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APPENDIX A 


TRANSPORT PROPERTIES OF AN AIR-HELIUM MIXTURE 

By using Wilke's approximations (ref. 21) the viscosity of a mixture of 1 
nents can be written as 


N 

^^mix " ^ ^i 


i=l 




N 


1 + 


2 

k=l ^ 

ym 


where 


<^ik " 




-l/2i 


1 



,J>h\ 

2^ 

“k 


' ► 

J 


“|2 


The expression for the thermal conductivity is expressed as 


^mix 


N 

= )^ki' 


i=l 


N 


1 + 


‘^ik 


k=l 


Ci 


where 


^;k= 



“ 

-1/2 

, .1/2/ a/4 

1.065 



ij^] 

2\l2 

^k 


W rv 
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For these computations, 


ks =4 9r= - 5 
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APPENDIX A 


The expression for the viscosity of helium is broken up into two temperature ranges. 
For T < 1000 K, 

juigg = 2.3479 X 10’ 

and for T g 1000 K, 

Mfle = 1-6215 X lO’^ + 1.02085 x 10~® T - 8.33348 x 10'^^ -h 7.5447 x 10'^^ 

which is a curve fit of the data of Amdur and Mason (ref. 22). 

For temperatures less than 3000 K, the viscosity of air is taken from Hansen 
7 rp3/2 

(ref. 23) as 8.1363 x 10" ' . For temperatures greater than 3400 K and less 

(T + 363.6) 

than 15 000 K, interpolation from Hansen's tables of ji = jm(T,p) is used to determine 
the viscosity. 
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APPENDIX B 


COMPUTED EIGENVALUES 


The computed eigenvalues for the case where heat conduction is considered are as 
follows: 

6\2Np^u2 + 4a2\2Np^ + 4XXNp^ - SNp^. + 3yXX i|/3(c2 - Ci) ^ C 2 + 

3Npr 2 2 


^1 = - 


6X 


^2 = - 


^NprU^ + 4a2x2Np^ + 4XXNp^ - 3Np^ + 3yXX i\|3(C2 - Cj) _ C 2 + 
3Npr ~ ^2 2 


Vn = 


6x2Np u2 + 4a2x% + 4XXN - 3Np^ + 3yXX 

. £_i — _j_ 1^2 + Cj 


3 3Npj. 

where X = 4/i/^Npj.Npg q 


r 


/ \ 

exp 

i 

f 8N|^X®a6 + 

18N|^x5xy + (24N^j. - 54Ny X^X - 288u2n|^X® 


+ <-27Np^x4x2y2 + 


144N^j. + 81Np IX^X^ - 216 U^n|j.X^X 


y + -48N^^ 


216N 


Pr, 


X^x2 + 


•288u2Np^ + 648 u2n|j. 


X^X) 


27x2x^y3 + 54Np^X^X^y2 


+ 72N^^x3x^y - 64 n|^ 
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/ 


^2 = i 

exp 

SN^p^xV . 

18N|^x5Xy + (24 n|^ - 54N|j.jy5x - 288 u2n^^X® 






+ <-27Npj.x4x2^ ^ 


1 


144n|^ + 8lNpjX%2 _ 216 u2n|^X^X 


y + -48 n;^^ 


- 216n|jx^X^ + (-288U%pj. + 648U^N|Jx®xJ>a2 - 27X^X^y3 + 54Np^X^X^y2 


+ 72N|^x3x3y - 64 n|^x3x3) /27n|j. 


+ Cg) /3 


and 


Cg = iX^ 4096a6x%p^U® + <j (-2048aSx® - 4096a%X® + 1024a4x2x'7Np^ + 


f 92 16a® 


3072a®y)xx® + (l3824a^ - 7680a^y)x^x4 


Npj. + f576a^y2 + 10368a4y 


- 15552a'^)x^X^N|^U^ + <||256a^®X® + 1024a®XX® + 512a®x2x'^ - 1024a^X®X®)Np^ 


{l68sfiy - 2304a®)xx® + (3840a®y - 8448a®)x^X^ + (3840a^y - 5376a^)x®X® 


+ (3072a2 - 1536a2y)x‘^x2 


iN|r^ 


288a®y2 - 1152a®y + 4320a®)x^x'^ + (2880a^y2 

(Equation continued on next page) 
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- 6336a4y + 10368a'^)x^x3 + ^4608a2y2 _ 3456a2y)x'^X^ + (-432a4y3 + 864a‘^y2 

- 3888a‘^y)x^x2 + (-864a2y3 - 2592 sl \^ x \^ Np^. + 1296a2y3x^X^U^ + (e4a®X^x‘^ 

+ 256a6x^X^ + 256a^X^x2)Npj, + ( 96 a® - 96a^y)x^X^ + (-192a®y - 192a6)x^X^ 

+ (384a^y - 1536a4)xV + (?68a2y - l536a2)x®X + (36a®y2 - 72a8y + 36aS)x^X^ 

+ (-144a6y2 + 720a®y - 576a6)x^X^ + ^-864aV + 2880a^y - 1152a^)xV + ^2304a2y 

- 576a2y2)x®X + 576y2x® + (l08a®y3 - 432a®y2 + 540a®y - 216a6jx^X^ 

+ (216a^y2 - 1080a‘^y2 + 864a4y)x‘^x2 + ^-432a2y3 - 432a2y2]x^X - 864y3x® Np^ 

\l/2 / 

+ (siaM - 162a4y3 + 81a4y^)x‘^x2 + (s24a2y4 - 324a2y3)x^X + 324y4x^ j / 3\/3 n|^ 
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TABLE 1.- EXPANSION TUNNEL INITIAL CONDITIONS 


Region 

p, mg/m3 

T, K 

U, m/s 

p, kN/m'2 

Gas 

0 

103.6449 

300 

0 

6.668 

Air 

0 

572.3 

3070.00 

2351.70 

512.3 

Air 

® 

9.7013 

1407.75 

5344.97 

3.0053 

Air 

@ 

0.0753 

300 

0 

0.02002 

Helium 


.1338 

300 

0 

.00860 

Air 

@) 

0.3927 

4935.00 

5344.97 

3.0053 

Helium 


1.838 

5373.30 

5344.97 

3.0053 

Air 


Uq = 5344.97 m/s Pq = 26.75 mg/m^ 

Rq = 886.47 N-sVk Xq = 0.0675 m 
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Figure 2.- Wave diagram, imperfect nozzle start. 
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Figure 3. 
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- Velocity and density distributions at succeeding time increments; 
= 16.67 N/m^; Mq = 5; no ternery diaphragm; inviscid case. 
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Figure 4.- Velocity and density distributions at succeeding time increments; = 16.67 N/m2; 

Mg ^ = 5; no ternery diaphragm; viscous case. 
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Figure 6.- Velocity, species concentration, and density distributions as a function of time; 
no ternery diaphragm; initial conditions given in table 1. 
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Figure 7.- x-t wave diagram of starting process for 10^ 

heJium acceleration gas. 






Time, microseconds 





1 nterface 

Secondary Starting Shock 

Downstream-Facing Characteristic 

Upstream-Facing Characteristic 


1 1 .1 1 1 1.1 

1 

X, meters 



Figure 10.- x-t wave diagram of starting process for 10° nozzle; 

= 0.1; helium acceleration gas. 
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Figure 11.- x-t wave diagram of starting process for 10® nozzle; 

0.01; helium acceleration gas. 






Figure 12.- x-t wave diagram of starting process for 8° nozzle; 

= 0*1; helium acceleration gas. 
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Figure 14.- x-t wave diagram of starting process for 10^ nozzle 
p^/p/CpN = 0.1; air acceleration gas. 
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Figure 15.- x-t wave diagram of starting process for lO^ozzle; 

= 0.01; air acceleration gas. 
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Figure 17.- Summary of available test time in 8° nozzle using helium acceleration gas 


as a function of A and 
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Figure 20,- Available test time as a function of inviscid core size and p 
in 10® nozzle for air and' helium acceleration gases. 
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